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Abstract 

Let s G (0, 1) be uniquely determined but only its approximations 
can be obtained with a finite computational effort. Assume one aims to 
simulate an event of probability s. Such settings are often encountered 
in statistical simulations. We consider two specific examples. First, the 
exact simulation of non- linear diffusions ([3]). Second, the celebrated 
Bernoulli factory problem ([10], [13]) of generating an f(p)— coin given a 
sequence X\, X2, ... of independent tosses of a p— coin (with known / and 
unknown p). We describe a general framework and provide algorithms 
where this kind of problems can be fitted and solved. The algorithms 
are straightforward to implement and thus allow for effective simulation 
of desired events of probability s. Our methodology links the simulation 
problem to existence and construction of unbiased estimators. 



1 Introduction 



Assume that one aims to simulate an event of unknown probability s £ (0,1) 
which is uniquely determined, however only its approximations can be obtained 
using a finite computational effort. Such settings are often encountered in sta- 
tistical simulations and emerge if e.g. s is given by a series expansion or a 
consistent estimator for s is available (see e.g. [7], [5], [5], [3] [TO], [12]). A 
celebrated example of this kind is the Bernoulli factory problem which moti- 
vated our work. It can be stated as follows. Let p 6 V C [0, 1] be unknown 
and let / : V — > [0, 1] . Then the problem is to generate Y, a single coin toss 
of an s — /(p)— coin, given a sequence X±,X2, ... of independent tosses of a 
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p— coin. For the historical context of this question and a range of theoretical 
results see [TB], [TO], [13], [H] and [PJ. In particular [T0J provide necessary and 
sufficient conditions for /, under which an algorithm generating an f(p)— coin 
exists. Nacu and Peres in [T3] suggest a constructive algorithm for simulating 
f(p) = min{2p, f — 2e} which is central to solving the problem for general / and 
allows for generating an f(p) — coin for a large class of functions (e.g. real ana- 
lytic, see [13] and Section [3] for details). The algorithm is based on polynomial 
envelopes of /. To run the algorithm one has to construct sets of {0, 1 } strings of 
appropriate cardinality based on coefficients of the polynomial envelopes. Un- 
fortunately its naive implementation requires dealing with sets of exponential 
size (we encountered e.g. 2 2 ) and thus is not very practical. Hence the authors 
provide a simple approximate algorithm for generating min{2p, 1}— coins. On- 
going research in Markov chain Monte Carlo and rejection sampling indicates 
that the Bernoulli factory problem is not only of theoretical interest (c.f. [JJ, [5], 
Chapter 16 of [2], and Section 0] of the present paper). However using approxi- 
mate algorithms in these applications perturbs simulations in a way difficult to 
quantify. 

In Section 2 we develop a framework for simulating events of unknown prob- 
abilities. Our approach is based on random sequences, say L n and U n under- 
and overestimating s that are monotone in expectations (i.e. E L n /* s and 
E [/„ \ s) and are reverse time super- and submartingales respectively. From 
L n and U n we construct L n and U n that are monotone almost surely and have 
the same expectations (c.f. Theorem 12. 51 Algorithm 4). Given L n and U n we 
sample events of probability s using a single U(0, 1) random variable. This result 
generalizes classical constructions for simulation of events of unknown probabil- 
ities using deterministic sequences (0). We link these results to existence and 
construction of unbiased estimators. In particular one can use the algorithms 
of Section [2] to obtain unbiased sequential estimators of a parameter of interest 
that is not necessarily in [0, 1]. 

We illustrate our results with examples. First, in Section [31 we present a 
reverse time martingale/unbiased estimator formulation of the Nacu-Peres al- 
gorithm which we believe gives a new perspective on the Bernoulli factory prob- 
lem. We identify the coefficients of the lower and upper polynomial envelopes 
as random variables of desired properties and implement the algorithm using 
a single U(0, 1) auxiliary random variable. We do not need to identify subsets 
of {0, 1} strings and thus avoid algorithmic difficulties of the original version. 
The martingale approach also simplifies the proof of validity of the Nacu-Peres 
algorithm. In the special case when / has an alternating series expansion with 
decreasing coefficients, the martingale approach results in new, very efficient 
algorithms. Second, in Section 2] we obtain the Exact Algorithm for diffusions 
introduced in [3] as an application of the generic Algorithm 3 of Section [2) 

2 Simulation of Events with Unknown Proba- 
bilities 

Throughout the paper we assume that we can generate uniformly distributed iid 
random variables Go, G%, ... ~ U(0, 1) which will serve as a source of randomness 
for algorithms. Thus to simulate an s— coin C s we just let C s := I{Gq < s}. We 
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will be concerned with settings where s is not known explicitly. 
The following simple observation will turn out very useful. 

Lemma 2.1. Sampling events of probability s £ [0, 1] is equivalent to construct- 
ing an unbiased estimator of s taking values in [0, 1] with probability 1. 

Proof. Let S, s.t. ES = s and F(S G [0, 1]) = 1 be the estimator. Then draw 
Go ~ C/(0, 1), obtain S and define a coin C s := 1{Gq < S}. Clearly 

P(C S = 1) = E I(G < S) = E (e (l(G < s) \ S = s)) = E5 = s. 

The converse is straightforward since an s— coin is an unbiased estimator of s 
with values in [0, 1]. □ 

Thus given S <E [0, 1], an unbiased estimator of s, we can sample events of 
probability s by the following algorithm. 

Algorithm 1. 

1. simulate Go ~ U(0, 1); 

2. obtain S; 

3. if Go < S set G s := 1, otherwise set G s := 0; 

4. output G s . 

Next assume that li,l2, ■■■ and Ui,it2, ... are sequences of lower and upper 
bounds for s converging to s. This setting is well known ([7j) and appears 
in a variety of situations, usually as an element of more complex simulation 
procedures, see e.g. [5], [15] . Here we use the following algorithm for simulating 
an s— coin. 

Algorithm 2. 

1. simulate Go ~ U(0, 1); set n = 1; 

2. compute l n and «„; 

3. if Go < l n set G s := 1; 

4. if G > u n set C s := 0; 

5. if l n < Go < u n set n := n + 1 and GOTO 2; 

6. output C s . 

The algorithm stops with probability 1 since l n and u n converge to s from 
below and from above. Precisely, the algorithm needs N > n iterations to stop 
with probability inf fc<„ Uk — sup fc< „ Ik- Because we can always obtain monotone 
bounds by setting u n '■= inffc<„ Uk and l n '■= sup fc< „ Ik, we assume that l n is an 
increasing sequence and u n is a decreasing sequence. 

The next step is to combine the above ideas and work with randomized 
bounds, i.e. in a setting where we have estimators L n and U n of the upper and 
lower bounds l n and u n - The estimators shall live on the same probability space 
and have the following properties that hold a.s. for every n = 1, 2, ... 

L n < U n (1) 
L n G [0, 1] and U n G [0, 1] (2) 
Ln-i < L n and U n -i > U n (3) 
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Note that we do not assume that L n < s or U n > s. Also condition © implies 
monotonicity of expectations, i.e. 

EL n = l n y S and E U n = u n \ s. (4) 

Let 

J"o = {0,fi}, T n = <r{L n , [/„}, J 7 k ,n=<j{Fk,Fk+i 1 ---Fn} for k < n. 
Consider the following algorithm. 
Algorithm 3. 

1. simulate Go ~ U(Q, 1); set n = 1; 

2. obtain L n and U n given To.n-i, 

3. if Go < L n set C s := 1; 

4. if Go > £4 set G s := 0; 

5. if L n < G < f/„ set n := ri + 1 and GOTO 2; 

6. output C s . 

Lemma 2.2. Assume (QP, (0), (Q) and |7p. Then Algorithm^ outputs a valid 
s—coin. Moreover the probability that it needs N > n iterations equals u n — l n . 

Proof. Probability that Algorithm [3] needs more then n iterations equals E(C/ n — 
L n ) — In — u n — > as n — > oo. And since < U n — L n is a decreasing sequence 
a.s., we also have U n — L n — > a.s. So there exists a random variable S 1 , such 
that for almost every realization of sequences {L n (u)} n >i and {U n (u))} n >i we 
have £„(w) / S(u) and [7 n (w) \ S(cj). By © we have S E [0, 1] a.s. Thus for 
a fixed cj the algorithm outputs an S(lu)— coin a.s. Clearly E L n < E S < E U n 
and hence ES=s. □ 

Remark 2.3. The random variable S constructed in the proof can be viewed as 
the unbiased estimator of s mentioned earlier with sequences L n and U n being 
its lower and upper random approximations. 

Remark 2.4. For Algorithm [3] assumption ([2|) can be relaxed to 



L n e (—oo,l] and U n G [0, 00) a.s. for every n = 1,2,... (5) 

The final step is to weaken condition ^ and let L n be a reverse time super- 
martingale and U n a reverse time submartingale with respect to !F n ,oo- Precisely, 
assume that for every n = 1, 2, ... we have 

E (L n -i I F n ,oo) = E (L n _i I J" n ) < L n a.s. and (6) 
E ([/„_! I ^„ :00 ) = E (Un-i I > U n a.s. (7) 



Consider the following algorithm, that uses auxiliary random sequences L 
and U n constructed online. 

Algorithm 4. 

1. simulate Go ~ U(0, 1); set n = 1; set i = ^0 = and Uq = Uq = 1 

2. obtain L„ and £/„ given J-o,n-i, 
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3. compute L* n = E (L n _i | T„) and [7 * = E | T n ). 

4. compute 



ir, 



(t/ n _i - L B _i) (9) 



f/* — r/ 

t/* - LI 



5. if Go < L n set C s := 1; 

6. if G > U n set G s := 0; 

7. if L n < G < U n set n := n + 1 and GOTO 2; 

8. output C s . 

Theorem 2.5. Assume (0), Q), @) and Then Algorithm^ outputs 
a valid s—coin. Moreover the probability that it needs N > n iterations equals 

Un ■ 

Proof. We show that L and U satisfy (P), ©, © and © and hence Algorithm 
2] is valid due to Lemma 12721 

Conditions (fTJ), © and © are straightforward due to construction of L and 
U and @, ©. 

To prove © we show that the construction in step 4 of Algorithm[4]preserves 
expectation, i.e. 

E L n = E L n = /„ and EU n — E U n = u n . (10) 

It is straightforward to check that (TIT))) holds for n = 1,2. Moreover note that 
C/o — io = 1 a.s., t/j* - L* = 1 a.s. and from © and © we have 

U n - L n = (Un-i - L n -i) — and hence 

^« - n - 1 + 77; — 7T77J r* 773 — j^{Ui-Li). (lij 

Now we compute E L n by induction, conditioning (|lip subsequently on ^"2,00, ■ ■ ■ , 
and using ([6]) and ©. Calculation of E U n is identical. 



ei„ = el^+eiei :;: ;::;r ;r x -^-^(^-^) 
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-1 
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T * TT* 
lJ n u n-l 
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-1 





Ui - L* 2 



2,00 



el, , - e 1 ;:- ' ; -' t— |e (U, - Ll \T 2 ^] 



= E Ln-i + E " " n - L n - L ■■■ (U 2 - L a ) 

" 1 \ TT* T* TT* T* TT* T* v z z/ 

= E L„_i + E (L n - L* n ) = E (E (L„-i + L n - £;|^„,oo)) = E L„. 

□ 

Remark 2.6. All of the discussed algorithms are valid if n takes values along an 
increasing sequence nj /* oo. 



5 



Now let us link once again the algorithmic development of this Section with 
construction of unbiased estimators. Lemma 12.11 together with Theorem 12.51 
result in the following construction of sequential unbiased estimators based on 
under- and overestimating reverse time super- and submartingale sequences. 
The estimators are sequential in the sense that the amount of input needed to 
produce them is random. 

Theorem 2.7. Suppose that for an unknown value of interest s £ R, there exist 
a constant M < oo and random sequences L n and U n s.t. 

P(L„ < U n ) = 1 for every n = 1, 2, ... 

P(L n £ [-M, M]) = 1 and V(U n £ [-M, M]) = 1 for every n = 1, 2, ... 

E L n = l n y s and E U n = u n \ s 

E (Ln-i | .Fn.oo) = E (i„_i | J"„) < L„ o.s. and 

E (£/„_! | jr„ >00 ) = E ([/„_! | > C/„ a.s. 

Then one can construct an unbiased estimator of s. 

Proof. After rescaling, one can use Algorithm 4 to sample events of probability 
(M+s) /2M, which gives an unbiased estimator of (M+s)/2M and consequently 
of s. □ 



3 Application to the Bernoulli Factory Problem 

Based on Section we provide here a practical version of the Nacu-Peres algo- 
rithm for simulating an f(j>)— coin from a sequence of p— coins, where f(p) = 
min{2p, 1 — 2e}. This is central to the general version of the Bernoulli factory 
problem, as [13] develops a calculus for collapsing simulation of a real analytic 
function, say g, to simulation of f(p) = min{2p, 1 — 2e}. Briefly, one takes a 
series expansion of g and uses a composition of appropriate techniques (e.g. for 
simulating a sum or a difference of simulable functions) . We refer to the original 
paper for details. 

In particular we prove Proposition 13 - 1 1 a general result, which is a minor 
modification of Proposition 3 in |13j . However its proof, different from the 
original one, links polynomial envelopes of / with the framework of Section [2] 
by identifying terms. It results in an immediate application of Algorithm [4j 

Proposition 3.1. An algorithm that simulates a function f onV C (0, 1) exists 
if and only if for all n > 1 there exist polynomials g n (p) o.nd h n (p) of the form 



— A- 



s.t. 

(i) < a(n,k) < b(n,k) < 1, 

(ii) lim^oo g n (p) = f{p) = lim^oo h n (p), 
(Hi) For all m < n, their coefficients satisfy 

k ( n ~ m ^ k ( TL ~ m ) ( m ) 

a(n,k) > k ~L\ 1 a(m,i), b(n,k) < ^ fc ~! 1 b{m,i). (12) 

i=0 \kJ i=0 \k) 
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Proof. We skip the implication algorithm => polynomials, as it has been shown 
in [13) , and focus on proving polynomials => algorithm using framework of Sec- 
tion O Let Xi, X 2 , ... be a sequence of independent tosses of a p— coin. Define 

random S6CJU6I1C6S {£/n? £^n}n 

>i as follows: if J27=i -^i = ^, tnen ^« — a ( n > ^) 
and U n = b(n, k). In the rest of the proof we check that |T|), (J5J), Q, (JH) and 
(JT]) hold for {L n ,U n } n >i with s = f(p). Thus executing Algorithm 2] with 
{L n , U n } n >\ yields a valid /(p)-coin. 

Clearly {TJ and © hold due to (i). For flU) note that E L„ = g n {p) / f(p) 
and EU„ = /i„(p) \ /(f>)- To obtain ([5]) and define the sequence of random 
variables H n to be the number of heads in {X\. .... X n J, i.e. H n = XaLi -^i 
and let £/„ = a{H n ). Thus L„ = a(n, H n ) and ?7 n = b(n,H n ), hence T n C (?„ 
and it is enough to check that E(L m \Q n ) < L n and E(Z7 m |C? n ) > U n for m < n. 
The distribution of iJ m given H n is hypergeometric and 

E(L m \g n ) = E(a(m, H m )\H n ) = ^ — jTrY~ a ( m ^) ^ a(n,H n ) = L n . 

i=o Iff J 

Clearly the distribution of -ff m given H n is the same as the distribution of H m 
given {H n , H n+ i, . . .}. The argument for U n is identical. □ 

Remark 3.2. In contrast to |13j . throughout this section we simulate / in the 
weak sense, i.e. we use U(0, 1) as an auxiliary random variable. This is the 
natural approach in applications and also this is equivalent to strong simulability 
if V C (0,1), c.f. [TO]. 

To give a more complete view of the Bernoulli factory problem in the frame- 
work of Section [21 and before moving on to practical aspects of the problem, we 
show, as a corollary from Lemma (23 a result originally established in [TU] and 
also provided in 13J, namely that generating min{2p, 1}— coins from p— coins is 
not possible. 

Corollary 3.3. An algorithm that simulates f(j>) — 2p for p G V = (0,1/2) 
does not exists. 

Proof. We show that there does not exists an unbiased estimator of 2p for p £ 
(0, 1/2) that takes values in [0, 1] and we conclude the corollary from Lemma l2~T1 
The idea of the proof is to show that such an estimator must take values smaller 
then 1/2 with strictly positive probability independent of p and then let p /* 1/2 
so that 2p / 1. 

Let S be such an estimator and let Xi, X 2 , . . . be a sequence of p— coins. We 
allow S to be sequential and use an auxiliary random variable Rq independent 
of the p— coins. So 

S = S({X±,X2, ■ ■ ■ ,},T,i? ) = S({Xi,X 2 , ■ ■ ■ ,Xt},Rq), 

where T is a stopping time with respect to cr{.?i ira , C/}, where {^ 7 i l n} n >i is the 
filtration generated by X\, X 2l ■ ■ ■ and Q is a cr— algebra independent of T\, n and 
generated by Rq. Clearly the joint distribution of {\X\ , X 2 , ■ ■ ■ }, T, Rq} depends 
on p. We denote it by V p and let P p u be the projection of P p on {Xi, X 2 , . . . , X t }. 
Now fix p = 1/4. Since 2p = 1/2 we have 5 := ¥ 1/4 {S < 1/2) > 0. Moreover 
there exists such an to that 

¥ 1/4 {S<l/2; T<t )>5/2. 
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Note that P p i to is absolutely continuous with respect to Pi/4|t for all p £ 
[1/4, 1/2) and 

inf inf ^^ttt > 2 _t °, 

pe[l/4,l/2)AC{0,l}*o P 1/4 | to (A) 

and consequently for every p £ [1/4, 1/2) we have 

F p (S < 1/2) > P p (5 < 1/2; T < t ) > 2- t °P 1/4 (S' < 1/2; T < t ) > 82- {U)+1) . 
Now let p y 1/2. This combined with S £ [0, 1] contradicts unbiasedness. □ 

Given a function /, finding polynomial envelopes satisfying properties re- 
quired by Proposition l3.1l is not easy. Section 3 of [13] provides explicit formulas 
for polynomial envelopes of f(p) = min{2p, 1 — 2e} that satisfy conditions of 
Proposition 13 - 1 1 precisely a(n, k) and b(n, k) satisfy (ii) and (iii) and one can 
easily compute uq = no(e) s.t. for n > no condition (i) also holds, which is 
enough for the algorithm (however no is substantial, e.g. no(e) = 32768 for 
e = 0.1 and it increases as e decreases). By Theorem 12.51 the probability that 
Algorithm [3] needs N > n inputs equals h n (p) — g n (p) . The polynomials pro- 
vided in (T3] satisfy h n {p) — g n (p) < Cp n for p £ [0,1/2 — 4e] guaranteeing 
fast convergence, and h n (p) — g n (p) < Dn^ 1 / 2 elsewhere. Using similar tech- 
niques one can establish polynomial envelopes s.t. h n (p) — g n (p) < Cp n for 
p £ [0, l]\(l/2 — (2 + c)e, 1/2 — (2 — c)e). We do not pursue this here, however 
in applications it will be often essential to obtain polynomial approximations 
tailored for a specific problem and with desired properties. Moreover, we note 
that despite the fact that the techniques developed in [T3] for simulating a real 
analytic g exhibit exponentially decaying tails, they arc often not practical. 
Nesting k times the algorithm for f(p) = min{2p, 1 — 2e} is very inefficient. 
One needs at least no(e) k of original p— coins for a single output. 

As mentioned earlier, a naive implementation of the Nacu-Peres algorithm 
requires dealing with sets of {0, 1} strings of exponential size. Other imple- 
mentations with reduced algorithmic cost are certainly possible with additional 
effort. However, our martingale approach that uses Algorithm 2] in the way indi- 
cated in the proof of Proposition [3T] avoids this problem completely (a C-code 
for f(p) = min{2p, 1 — 2e} is available on request). 

Nevertheless, we note that for both algorithms, i.e. the original Nacu-Peres 
version and our martingale modification, the same number of original p— coins 
will be used for a single f(p)— coin output with f(p) — min{2p, 1 — 2e} and 
consequently also for simulating any real analytic function using methodology 
of [13] Section 4. A significant improvement in terms of p— coins can be achieved 
only if the monotone super/sub-martingales can be constructed directly and 
used along with Algorithm [3] This is discussed in the next subsection. 

3.1 Bernoulli Factory for alternating series expansions 

In the following Proposition we describe an important class of functions for 
which an f(p)— coin can be simulated by direct application of Algorithm 3. 

Proposition 3.4. Let f : [0, 1] — » [0, 1] have an alternating series expansion 

oo 

f(p) = J2(-l) k a k p k with 1 > ffl > ai > • • ■ 

k=0 
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Then an f(p) — coin can be simulated by Algorithm 3 and the probability that it 
needs N > n iterations equals a n p n . 

Proof. Let Xi, X2, ... be a sequence of p— coins and define 

U ■= a L a := 0, 

U n -i - a n rifc=i x k if n is odd, 
L„_i if n is even, 



U n -i if n is odd, 

+ a n rifc=i -Xfc if n is even. 



Clearly |T|), J2]), © and J4} are satisfied with s = f(p). Moreover, 
u n - l n = 1 C/i, - IE in = a«P™ < On- 

Thus if a n — > 0, the algorithm converges for p 6 [0, 1], otherwise for p £ [0, 1). 

□ 

Next we illustrate the difference between application of Algorithm 2] based 
on the Nacu-Peres approach and direct usage of Algorithm [3] for simulating 
f(p) — exp(— ap), a < 1. This function appears in applications discussed in 
Section [5] Weak simulation is considered (c.f. Remark 13.21 and [TU]). i.e. all 
normally available random variables are obtained directly, not from p— coins. 

First consider sampling an f(p)— coin by collapsing the problem to doubling 
(i.e. sampling of min{2p, 1 — 2s}) using techniques of [13 Section 4 and Algo- 
rithm[H we refer to the original paper for a complete description of the approach. 
The aim is to use as few doubling steps as possible, since doubling is expensive. 
Let k :=€ {1, 2} be s.t. 2 k > e a . We have 

00 n 2n 00 „2n+l 

-ap = \ ^ _^ 2n _ \ ^ " „2n+l 

^ (2n)r ^ (2n+l)V 

n=0 v ' n=0 v ; 
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rj(*+(p)-*-(p))\ as) 

e -a a 2n ~ e -a a 2n+l 

where . + (p) = and '-to = E^+ijT^- 

First consider obtaining (s + (p) — s_ (p)) —coins. This will be done by reversing 
(l — s + (pj) + s_(p)j— coins. Since 



1 



-(p) + S _(p) = 2(i(l- a+ (p)) + i«_(p)), (14) 



one feeds the doubling algorithm with (|(1 — s+(p)) + ^s_(p)) — coins obtained 
by tossing a fair coin first and using an (1 — s+(p)) — coin or an s_(p) — coin in 
the second step, depending on the outcome of the fair coin. 

We now describe sampling an s + (p)— coin C s+ ( p ). An s_(p)— coin can be 
obtained in a similar manner. Due to the specific form of series expansion 
using a Poisson mixture is more efficient then using the enforced geometric 
mixture suggested in the proof of Proposition 16 [13j . details are below. Sample 
N ~ Poiss(a). If N is even then generate iid p— coins X\, . . . , Xn and declare 
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C s+ ( p ) := 1 if Xi = ■ ■ ■ = Xn = 1- Otherwise, i.e. if N is odd or N is even 
and 3i<k<N s.t. Xf. = 0, declare C s+ ^ := 0. Finally, ^—thinning should be 
applied to the (s + (p) — s_(p))— coins and the doubling algorithm should be 
nested k times on |r(s + (p) — s„(p)) —coins. 

To approximate the total simulation effort let K p be the cost of obtaining 
the p— coin and assume U(0, 1) r.v's cost to be 0(1). Then the cost of a fair 
coin and the Poisson r.v. is also O(l) (c.f. [7]). We assume K p ^> 1. Since 
s+(p) - s_(p) = e-^P+i) > e - 2a , 

\(1- S+ (p)) + ^-(p) < \~\ e ~ 2a ' 

and in (fT4|) one can take e\ = e _2a /5 for the doubling algorithm for min{2p, 1 — 
2ei}. Moreover to apply k times the doubling scheme in (fT3"|) we have to ensure 
that 8 < e~ ap < 1 — 6 (c.f. Proposition 17 of [T3]) and therefore we have 
to restrict our considerations to the situation where a lower bound on p, say 
pi is known. This implies that f-(s + (p) - s_(p)) = e" ap /2 < e"^' /2 = 
h — (h — e~ api /2) and in the last iteration of the doubling scheme in (fT3|) one 
can take say £2 = (h — e~ api /2)/3. This yields a lower bound on the number of 
p— coins required before the algorithm can stop (c.f. the discussion in Section[3]), 
namely n.o(^i)"-o(£2)- If fc = 2, the bound must be multiplied by no(e) with e 
used in the first iteration of the doubling scheme in (fT5)l . Recall that e.g. 
no(e) — 32768 for e = 0.1 and it increases as e decreases. Therefore, when 
applying Algorithm 4 based on the Nacu-Peres approach, a conservative lower 
bound on the total simulation effort is 2 30 K p . 

On the other hand, for the exponential function we readily have an alter- 
nating series expansion and can apply Algorithm 3 directly by appealing to 
Proposition [23 Then, we have 

F{N > n) = Ml, 

TV. 

where N is the number of p— coins required. This implies E7V < e, and the 
expected simulation effort is bounded from above by eK p and the bound holds 
uniformly for a G [0, 1] and p E [0, 1]. 

4 Application to the exact simulation of diffu- 
sions 

In this Section we derive the Exact Algorithm for diffusions introduced by [3] 
as a specific application of the Bernoulli factory for alternating series of Section 
13.11 We are interested in simulating Xt which is the solution at time T > of 
the following Stochastic Differential Equation (SDE): 

dX t = a{X t ) dt + dWt, Xa = xeR,te [0, T] (15) 

driven by the Brownian motion {Wt ; < t < T}, where the drift function a 
is assumed to satisfy the regularity conditions that guarantee the existence of 
a weakly unique, global solution of (fT5"|) . see ch.4 of [TT]. Let fl = C([0,T],R) 
be the set of continuous mappings from [0, T] to R and lu be a typical element 
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of f2. Consider the co-ordinate mappings B$ : O i— *• R, t £ [0, T], such that for 
any t, B t {uj) — uj(t) and the cylinder er-algebra C — cr({B t ; < t < T}). We 
denote by W x — {W* ; < t < T} the Brownian motion started at x £ R, and 
by W x ' u — {Wt' u ; < t < T} the Brownian motion started at x and finishing 
at u £ R at time T; the latter is known as the Brownian bridge. We make the 
following assumptions for a: 

1. The drift function a is differentiable. 

2. The function h(u) = cxp{A(u) - (u - x) 2 /2T}, u £ R, for A(u) = 
Jo a {y)^Vi is integrable. 

3. The function (a 2 + a )/2 is bounded below by I > -co, and above by 
r + £ < oo. 

Then, let us define 

4>( u ) = -[(a 2 + a')/2 -I] £ [0, 1] , (16) 
r 

Q be the probability measure induced by the solution X of (I15|) on (f2,C), W 
the corresponding probability measure for W x , and 1 be the probability mea- 
sure defined as the following simple change of measure from W: dW/dZ(w) oc 
cxp{— A(Bt)}- Note that a stochastic process distributed according to Z has 
similar dynamics to the Brownian motion, with the exception of the distribution 
of the marginal distribution at time T which is biased according to A. Hence, 
we refer to this process as the biased Brownian motion. In particular, the biased 
Brownian motion conditional on its value at time T has the same law as the 
corresponding Brownian bridge. 

The final steps of the mathematical developement entail resorting to the 
Girsanov transformation of measures (see for instance ch.8 of [T3]) to obtain 
dQ/dW; applying an integration-by-parts (possible by means of Assumption 1) 
to eliminate the stochastic integral involved in the Radon-Nikodym derivative; 
and using the definition of Z to obtain that 

^ (w) oc exp | - rT T' 1 <j>(B t )dt^ < 1 Z - a.s. (17) 

The details of this argument can be found in [31 [5]. By a standard rejection sam- 
pling principle, it follows that a path uj generated according to Z and accepted 
with probability (fl7|) . yields a draw from Q. Hence, the following algorithm 
yields an exact sample from the solution of (| 1 5[) at time T: 

1. simulate u ~ h 

2. generate a C s coin where s :— e~ rTJ , and J := f Q T~ 1 (f>(Wt' u )dt; 

3. If C s = 1 output u and STOP; 

4. If C s = GOTO 1. 

Exploiting the Markov property, we can assume from now on that rT < 1. If 
T is such that rT > 1, then we can devise sub-intervals of length 8 such that 
rS < 1 and apply the algorithm sequentially. 
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Clearly, the challenging part of the algorithm is Step 2, since exact compu- 
tation of J is impossible due to the integration over a Brownian bridge path. 
On the other hand, it is easy to generate J-coins: Cj — !(?/> < (j)(W£' u )), where 
tp ~ U(0, 1) and x ~ U(0,T) independent of the Brownian bridge W x ' u and of 
each other. Therefore, we deal with another instance of the problem studied 
in this article: given p-coins how to generate /(p)-coins, where here / is the 
exponential function. This is precisely the context of Section 13.11 where the 
use of Algorithm [3] was advocated for efficient simulation. As a final remark, 
we note that exact simulation algorithms have been proposed in [H [5] for mul- 
tivariate diffusions and unbounded drit functionals. These extensions involve 
decompositions of the Brownian motion and auxiliary Poisson processes. 
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